Random Developement of  Regular hexadecachoron
(To see how the program work, please refer to “working”.)

Animate[

OA={1/Sqrt[2],-1/Sqrt[6],-Sqrt[3]/6};
OB={0,Sqrt[6]/3,-Sqrt[3]/6};
OC={-1/Sqrt[2],-1/Sqrt[6],-Sqrt[3]/6};
OD={0,0,Sqrt[3]/2};
dpts=cpts=pts0={N[{OA,OB,OC,OD}]};


(*get the next polyhedron; For example n=4, x={A,B,C,D} then OE=2/3(OA+OB+OC)-OD*)
getNext[x_List,n_Integer]:=ReplacePart[x,n->(2/3 Total[Drop[x,{n}]]-x[[n]])];

(*cpts=continuous developement's 3Dpoints set*)
For[i=1,i≤15,i++,
AppendTo[cpts,getNext[cpts[[i]],nums[[i]]]]
];

(*dpts=discontinuous developement's 3Dpoints set*)
For[i=1,i≤15,i++,
If[
Count[mycells[[i]],2]≠1,AppendTo[dpts,getNext[dpts[[i]],nums[[i]]]],
(m=Position[mycells[[i]],2][[1,1]];
AppendTo[dpts,getNext[pts0[[1]],m] ])
]
];

(*"snake" is a continuous developement like a snake. "starfish" is a discontinuous developement.*)

Switch[type,"snake",points=cpts,"starfish",points=dpts];

(*x,y,z range*)
Xcords=Flatten[points,1][[All,1]];
Ycords=Flatten[points,1][[All,2]];
Zcords=Flatten[points,1][[All,3]];
(*{Xcords,Ycords,Zcords}=Table[Table[Table[points[[j,k,cord]],{k,1,4}],{j,1,16}],{cord,1,3}];*)
{Mx,mx}={Max[Xcords],Min[Xcords]};
{My,my}={Max[Ycords],Min[Ycords]};
{Mz,mz}={Max[Zcords],Min[Zcords]};


(*definition of tetrahedron*)
tetrahedron[vecs_]:=Table[Polygon[Drop[vecs,{k}]],{k,1,4}];

(*contents of Graphics3D*)
mytext=If[showname,Text[Style[StringJoin["No",ToString[Floor[t]+1],". ", names[[Floor[t]+1]]],15],{mx,my,mz},{-1,-1}],Point[{100,0,0}]];
polyABCD={Opacity[0.9],Red,tetrahedron[points[[1]]]};
polys=Table[If[k≠Floor[t]+1,{Opacity[0.9],tetrahedron[points[[k]]]},
{Opacity[0.9(t-Floor[t])],tetrahedron[points[[k]]]}],{k,2,Floor[t]+1}];

(*show*)
Graphics3D[{polyABCD,polys,mytext},PlotRange->{{mx,Mx},{my,My},{mz,Mz}},Boxed->False
]
,
(*Animation Control*)
{t,0,16-10^(-10)},{type,{"starfish","snake"},ControlType->SetterBar},
{{showname,True,"show name"},{True,False},ControlType->Checkbox},

(*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*)
(*Very loooong intialization. To set up "nums" & "mycells", RANDOMLY.*)
Initialization:>(cells=cells0={(*{1,1,1,2},*){1,1,2,1},{1,1,2,2},{1,2,1,1},{1,2,1,2},{1,2,2,1},{1,2,2,2},{2,1,1,1},{2,1,1,2},{2,1,2,1},{2,1,2,2},{2,2,1,1},{2,2,1,2},{2,2,2,1},{2,2,2,2}};

mycells0={{1,1,1,2}};

(*create "mycell" and "num". "mycell" is aligned such that each member is different from the next member only by one part, which is shown by "num". *)

While[ Length[cells]≥1,
Catch[(
cells=cells0;
mycells=mycells0;
nums={4};(*mycells starts from {1,1,1,2}. And the base cell is {1,1,1,1}*)
For[k=1,k≤15,k++,
(
i=1;
spl=Append[RandomSample[Range[4],4],0];(*for example, spl={2,3,1,4}*)
n=spl[[1]];(*for example, n=2*)
While[(j=Position[cells,nextcell=ReplacePart[mycells[[k]],n->3-mycells[[k,n]] ]])=={},(i++;n=spl[[i]];If[i==5,Throw[sippai]] ) ];(*the main part of this func*)
cells=Delete[cells,j];
AppendTo[mycells,nextcell];
AppendTo[nums,n];
)](*end of for*)
)](*end of catch*)
];

(*partitioned mycells. Partitined by {1,1,1,2},{1,1,2,1},{1,2,1,1},{2,1,1,1}*)
parcells=(threeOne=Append[Position[Table[Count[mycells[[k]],2],{k,1,15}],1],{16}];Table[Take[mycells,{threeOne[[k,1]],threeOne[[k+1,1]]-1}],{k,1,4}]);

getname[x_]:=StringJoin[{x[[1]]/.{1->"A",2->"a"},x[[2]]/.{1->"B",2->"b"},x[[3]]/.{1->"C",2->"c"},x[[4]]/.{1->"D",2->"d"}}];

names=Prepend[ Map[getname,mycells],"ABCD"];


)
]

parcells//MatrixForm
nums